---
title: "f3 and life-history traits"
author: "Pierre Barry"
date: "`r format(Sys.time(), '%d %B, %Y, %H:%M')`"
output:
flexdashboard::flex_dashboard:
theme: paper
orientation: rows
social: menu
source_code: embed
vertical_layout: scroll
---
```{r global, include=FALSE}
list.of.packages <- c("readxl","ggplot2","ggpubr","ggrepel","plotly","My.stepwise")
for (i in list.of.packages){
if (i %in% installed.packages()[,"Package"] == FALSE){
install.packages(i)
}
eval(bquote(library(.(i))))
}
```
f3 and life-history traits {data-icon="ion-ios-pulse"}
=======================================================================
```{r}
data = read.table(paste(path,"/output/pop_genomics/f3.csv",sep=""),header = TRUE, sep=",",dec=".")
data = data[,-1]
lfh = as.data.frame(read_excel(paste(path,"/data/lfh/lfh.xlsx",sep="")),sheet=1)
lfh = lfh[,-1]
colnames(lfh)[2] = "Species"
data = merge(lfh,data,on=c("Species"))
fst = read.table(paste(path,"/output/pop_genomics/fst.csv",sep=""),header = T, sep=",")
fst = fst[,-1]
fst = fst[order(fst$Species), ]
data = merge(data,fst,on=c("Species"))
data$italic_species = italic_species=c("italic('A. boyeri')",
"italic('A. fallax')",
"italic('C. galerita')",
"italic('C. julis')",
"italic('D. labrax')",
"italic('D. puntazzo')",
"italic('G. niger')",
"italic('H. guttulatus')",
"italic('L. budegassa')",
"italic('L. mormyrus')",
"italic('M. merluccius')",
"italic('M. surmuletus')",
"italic('P. erythrinus')",
"italic('S. cabrilla')",
"italic('S. cantharus')",
"italic('S. cinereus')",
"italic('S. pilchardus')",
"italic('S. sarda')",
"italic('S. typhle')")
admixture_f3_significant = c()
admixture_f3_negative = c()
admixture_f3_negative_number = c()
for (i in 1:nrow(data)){
for (j in 1:ncol(data)){
if (is.na(data[i,j])==TRUE){
data[i,j] = 0
}
}
}
for (i in 1:nrow(data)){
if (data$f3_MU_LI_GA[i] < 0 | data$f3_FA_LI_GA[i] < 0){
admixture_f3_negative[i]="Yes"
admixture_f3_negative_number[i] = 1
} else {
admixture_f3_negative[i]="No"
admixture_f3_negative_number[i] = 0
}
if (data$f3_MU_LI_GA_Z[i] < -3 | data$f3_FA_LI_GA_Z[i] < -3){
admixture_f3_significant[i]="Yes"
} else {
admixture_f3_significant[i]="No"
}
}
data$admixture_f3_negative = admixture_f3_negative
data$admixture_f3_negative_number = admixture_f3_negative_number
data$admixture_f3_significant = admixture_f3_significant
```
Row {.tabset data-height=1000}
-----------------------------------------------------------------------
```{r}
plot_legend =c("Body size (cm)",
"Trophic level",
"log[Fecundity (eggs/day)]",
"Propagule size (mm)",
"Age at first maturity (years)",
"Lifespan (years)",
"Adult lifespan (years)",
"Pelagic Larval Duration (days)",
"Hudson's Fst - Gulf of Lion - Costa Calida",
"Hudson's Fst - Gulf of Lion - Algarve",
"Hudson's Fst - Gulf of Lion - Bay of Biscay",
"Hudson's Fst - Costa Calida - Algarve",
"Hudson's Fst - Costa Calida - Bay of Biscay",
"Hudson's Fst - Algarve - Bay of Biscay"
)
my_comparisons <- list( c("No","Yes"))
```
### Body Size
```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
x1=data$admixture_f3_negative_number,
y=data$Body_Size,
Species = data$Species_plot,
italic_species = italic_species)
data_tmp = data_tmp[data_tmp$y!="NA",]
model<-glm(x1 ~ y,
data = data_tmp,
family=binomial)
p<-ggplot(data_tmp,aes(x = x, y = y)) +
geom_violin(alpha = 0.1,aes(fill=x),col="white") +
geom_jitter(aes(text=paste("Admixture: ", x, "\n Body size (cm) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
theme_classic() +
#labs_pubr() +
xlab("Presence of admixture indicated by negative f3") +
ylab("Body size (cm)") +
geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
label = "p.format") +
scale_fill_viridis_d(begin = 0.25,end=0.75) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
coord_flip()
ppp<-ggplotly(p,height=500,tooltip=c("text")) %>% partial_bundle()
ppp
```
### Trophic level
```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
x1=data$admixture_f3_negative_number,
y=data$Trophic_Level,
Species = data$Species_plot,
italic_species = italic_species)
data_tmp = data_tmp[data_tmp$y!="NA",]
model<-glm(x1 ~ y,
data = data_tmp,
family=binomial)
p<-ggplot(data_tmp,aes(x = x, y = y)) +
geom_violin(alpha = 0.1,aes(fill=x),col="white") +
geom_jitter(aes(text=paste("Admixture: ", x, "\n Trophic level : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
theme_classic() +
#labs_pubr() +
xlab("Presence of admixture indicated by negative f3") +
ylab("Trophic level") +
geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
label = "p.format") +
scale_fill_viridis_d(begin = 0.25,end=0.75) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
coord_flip()
ppp<-ggplotly(p,height=500,tooltip=c("text")) %>% partial_bundle()
ppp
```
### Propagule size
```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
x1=data$admixture_f3_negative_number,
y=data$Propagule_Size,
Species = data$Species_plot,
italic_species = italic_species)
data_tmp = data_tmp[data_tmp$y!="NA",]
model<-glm(x1 ~ y,
data = data_tmp,
family=binomial)
p<-ggplot(data_tmp,aes(x = x, y = y)) +
geom_violin(alpha = 0.1,aes(fill=x),col="white") +
geom_jitter(aes(text=paste("Admixture: ", x, "\n Propagule size (cm) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
theme_classic() +
#labs_pubr() +
xlab("Presence of admixture indicated by negative f3") +
ylab("Propagule size (mm)") +
geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
label = "p.format") +
scale_fill_viridis_d(begin = 0.25,end=0.75) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
coord_flip()
ppp<-ggplotly(p,height=500,tooltip=c("text")) %>% partial_bundle()
ppp
```
### Fecundity
```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
x1=data$admixture_f3_negative_number,
y=data$Fec,
Species = data$Species_plot,
italic_species = italic_species)
data_tmp = data_tmp[data_tmp$y!="NA",]
data_tmp$y = log(as.numeric(data_tmp$y))
model<-glm(x1 ~ y,
data = data_tmp,
family=binomial)
p<-ggplot(data_tmp,aes(x = x, y = y)) +
geom_violin(alpha = 0.1,aes(fill=x),col="white") +
geom_jitter(aes(text=paste("Admixture: ", x, "\n Fecundity (eggs/day): ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
theme_classic() +
#labs_pubr() +
xlab("Presence of admixture indicated by negative f3") +
ylab("log[Fecundity (eggs/day)]") +
geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
label = "p.format") +
scale_fill_viridis_d(begin = 0.25,end=0.75) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
coord_flip()
ppp<-ggplotly(p,height=500,tooltip=c("text")) %>% partial_bundle()
ppp
```
### Age at first maturity
```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
x1=data$admixture_f3_negative_number,
y=data$Age_Mat,
Species = data$Species_plot,
italic_species = italic_species)
data_tmp = data_tmp[data_tmp$y!="NA",]
model<-glm(x1 ~ y,
data = data_tmp,
family=binomial)
p<-ggplot(data_tmp,aes(x = x, y = y)) +
geom_violin(alpha = 0.1,aes(fill=x),col="white") +
geom_jitter(aes(text=paste("Admixture: ", x, "\n Age at first maturity (years) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
theme_classic() +
#labs_pubr() +
xlab("Presence of admixture indicated by negative f3") +
ylab("Age at first maturity (years)") +
geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
label = "p.format") +
scale_fill_viridis_d(begin = 0.25,end=0.75) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
coord_flip()
ppp<-ggplotly(p,height=500,tooltip=c("text")) %>% partial_bundle()
ppp
```
### Lifespan
```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
x1=data$admixture_f3_negative_number,
y=data$Lifespan,
Species = data$Species_plot,
italic_species = italic_species)
data_tmp = data_tmp[data_tmp$y!="NA",]
model<-glm(x1 ~ y,
data = data_tmp,
family=binomial)
p<-ggplot(data_tmp,aes(x = x, y = y)) +
geom_violin(alpha = 0.1,aes(fill=x),col="white") +
geom_jitter(aes(text=paste("Admixture: ", x, "\n Lifespan (years): ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
theme_classic() +
#labs_pubr() +
xlab("Presence of admixture indicated by negative f3") +
ylab("Lifespan (years)") +
geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
label = "p.format") +
scale_fill_viridis_d(begin = 0.25,end=0.75) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
coord_flip()
ppp<-ggplotly(p,height=500,tooltip=c("text")) %>% partial_bundle()
ppp
```
### Adult lifespan
```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
x1=data$admixture_f3_negative_number,
y=data$Adult_Lifespan,
Species = data$Species_plot,
italic_species = italic_species)
data_tmp = data_tmp[data_tmp$y!="NA",]
model<-glm(x1 ~ y,
data = data_tmp,
family=binomial)
p<-ggplot(data_tmp,aes(x = x, y = y)) +
geom_violin(alpha = 0.1,aes(fill=x),col="white") +
geom_jitter(aes(text=paste("Admixture: ", x, "\n Adult lifespan (years) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
theme_classic() +
#labs_pubr() +
xlab("Presence of admixture indicated by negative f3") +
ylab("Adult lifespan (years)") +
geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
label = "p.format") +
scale_fill_viridis_d(begin = 0.25,end=0.75) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
coord_flip()
ppp<-ggplotly(p,height=500,tooltip=c("text")) %>% partial_bundle()
ppp
```
### Pelagic larval duration
```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
x1=data$admixture_f3_negative_number,
y=data$PLD,
Species = data$Species_plot,
italic_species = italic_species)
data_tmp = data_tmp[data_tmp$y!="NA",]
model<-glm(x1 ~ y,
data = data_tmp,
family=binomial)
p<-ggplot(data_tmp,aes(x = x, y = y)) +
geom_violin(alpha = 0.1,aes(fill=x),col="white") +
geom_jitter(aes(text=paste("Admixture: ", x, "\n Pelagic larval duration (days) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
theme_classic() +
#labs_pubr() +
xlab("Presence of admixture indicated by negative f3") +
ylab("Pelagic Larval Duration (days)") +
geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
label = "p.format") +
scale_fill_viridis_d(begin = 0.25,end=0.75) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
coord_flip()
ppp<-ggplotly(p,height=500,tooltip=c("text")) %>% partial_bundle()
ppp
```
### Model selection
```{r}
My.stepwise.glm(Y = "admixture_f3_negative_number",
variable.list = c("Body_Size","Trophic_Level","Propagule_Size",
"Age_Mat","Lifespan","Adult_Lifespan",
"Hermaphrodism","Parental_Care","PLD"
#,"FST_HUDSON_LI_MU","FST_HUDSON_LI_FA","FST_HUDSON_LI_GA",
#"FST_HUDSON_MU_FA","FST_HUDSON_MU_GA","FST_HUDSON_FA_GA"
),
data = data, myfamily = "binomial")
```
Packages used {data-icon="fa-map"}
=======================================================================
```{r}
installed.packages()[names(sessionInfo()$otherPkgs), "Version"]
```